librsb

A Shared Memory Parallel Sparse Matrix computations Library for the Recursive Sparse Blocks Format


Introduction

librsb is a library for sparse matrix computations implementing the (RSB) matrix format for multicore-parallel, shared-memory systems, via OpenMP.
It provides the most common operations for iterative solvers, e.g.: The RSB format is especially well suited for symmetric matrices (A ⩵ AT) and transposed (y ← β y + α AT B) multiplication variants.
librsb also implements the Sparse BLAS standard, as specified in the [BLAS Technical Forum] documents.

You might be interested in:

The Recursive Sparse Blocks (RSB) Format

The RSB data structure is hierarchical and allows cache-efficient and multi-threaded (that is, shared memory-parallel) operations on large sparse matrices. Its basic idea is a recursive subdivision of a matrix in the two dimensions. The subdivision algorithm attempts to partition the matrix until the individual submatrices occupy approximately the amount of memory contained in the CPU caches. The number of cores/thread available also plays a role.
The following picture shows the symmetric matrix audikw_1 (943695 rows, 39297771 nonzeroes, obtained from the SuiteSparse Matrix Collection and partitioned for different numbers of executing threads (respectively 1,2,4 and 8). Because of the symmetry, the upper triangle representation can be omitted.

Example of matrix partitioning for 1 thread.

Example of matrix partitioning for 2 threads.

Example of matrix partitioning for 4 threads.

Example of matrix partitioning for 8 threads.

The memory layout of the individual sparse submatrices follows a Z (or Morton) ordering. Each submatrix is stored either in coordinate or CSR (Compressed Sparse Rows) format, and the numerical index type size is chosen in order to minimize memory consumption.
This layout enables efficient and multithreaded algorithms implementing the Sparse BLAS operations.
The algorithms and the performance of RSB in practice are described in a series of papers; give a look at the papers .

Download

The following distributions or projects package librsb, according to their own specific method:

User Documentation & Literature

Reference Documentation and README

The complete documentation in HTML format for the latest released 1.3 version is available here.
The same documentation is included in the sources archive and provides several examples for starting using librsb.
For the old version 1.2, click here.

Papers, Posters, Presentations


Performance Guidelines

librsb is efficient on large matrices; that is, matrices whose memory occupation exceeds the memory cache. On current machines, this boundary is around a few million of nonzeroes. You might find librsb slower than a good CSR/CSC serial implementation if your matrices are smaller than that; but also faster if larger than that.
Keep in mind that what contributes the RSB format in librsb to be fast is: symmetric matrices, large matrices, many right hand sides, multiple cores.
You are welcome to give feedback about any performance related doubt, or provide a test case.

Build Instructions (tested on GNU/Linux)

Building the librsb library, C headers, documentation

Most numerical kernels code is auto-generated, and the supported numerical types can be chosen by the user at build time.
librsb can also be built serially (without OpenMP parallelism), if required.
So there are many options to build a customized library (see documentation).
However, building with reasonable defaults can be as simple as: 
# unpack the sources archive first; then:
./configure --help # there are many options available
./configure --enable-fortran-module-install --enable-matrix-types="double, double complex" \
	CC=gcc CXX=g++ FC=gfortran CFLAGS=-O3 CXXFLAGS=-O3 FCFLAGS=-O3 --prefix=$HOME/local
# If using other compilers, or a combination of different ones, you may have to use custom LD=..
# With GCC, you may get faster code with "-O3 -march=native -mfpmath=sse -msse2"...

make cleanall      # cleans up generated code
make               # builds the library and examples
make qtests        # quick consistency test
make tests         # long consistency test
make install       # install to $HOME/local
Please note that if you intend to link librsb to shared libraries (e.g.: in case of sparsersb's sparsersb.oct) then in the above you should specify flags to build position indipendent code, e.g.: CFLAGS="-O3 -fPIC" FCFLAGS="-O3 -fPIC" .
If you want to have a debug build with error verbosity activated, you may configure with something like: 
# good for familiarizing and debugging:
./configure --enable-allocator-wrapper --enable-interface-error-verbosity="2" \
 --enable-internals-error-verbosity="1" CFLAGS="-O0 -ggdb" CXXFLAGS="-O0 -ggdb" FCFLAGS="-O0 -ggdb" #...

Building the sparsersb Package for GNU Octave

 Assume you have both of librsb and a recent GNU Octave (complete with development programs) installed. Then you can build sparsersb by: 
% Be sure you compiled librsb with the position independent code flag (-fPIC).
% Browse to the sparsersb web page and get the sources archive, e.g.: sparsersb-1.0.9.tar.gz
% Start Octave:
% $ octave
% From Octave's prompt:
more off; pkg -local -verbose install sparsersb-1.0.9.tar.gz
% ... 
pkg load sparsersb
% ... 
test sparsersb
% ... 
help sparsersb
% ...
Alternatively, in the case you don't have a librsb installation at hand, you can have sparsersb use a librsb sources archive and build all in one go (linking statically): 
% First, download and save a librsb archive, say in ~/librsb-1.3.0.1.tar.gz
% Start Octave with environment variable LIBRSB_TARBALL set to the archive location:
% LIBRSB_TARBALL=$HOME/librsb-1.3.0.1.tar.gz octave
% or use setenv from within Octave:
setenv "LIBRSB_TARBALL" ~/librsb-1.3.0.0.tar.gz

% then from Octave's prompt:
more off; pkg -local -verbose install sparsersb-1.0.9.tar.gz
% this will take some more time
% ... 
pkg load sparsersb
% ... 
test sparsersb
% ... 
help sparsersb
% ...
Download instructions for both packages are here.
If you think any information here is wrong or out of date, consider reporting; see the contact section for this.

Usage Examples

The sparsersb Package for GNU Octave

The sparsersb package for GNU Octave allows you to manipulate sparse matrices using librsb without the need of programming in C/C++ or Fortran.
You can reuse your GNU Octave programs employing sparse matrices, and instead of creating matrices using the "sparse" keyword, use "sparsersb", like in this example: 
% from the Octave prompt:

R=(rand(3)>.6)
%  R =
%  
%     0   0   0
%     0   0   0
%     1   0   1

A_octave=sparse(R)
% A_octave =
% 
% Compressed Column Sparse (rows = 3, cols = 3, nnz = 2 [22%])
% 
%   (3, 1) ->  1
%   (3, 3) ->  1

A_librsb=sparsersb(R)
% A_librsb =
% 
% Recursive Sparse Blocks  (rows = 3, cols = 3, nnz = 2 [22%])
% 
%   (3, 1) -> 1
%   (3, 3) -> 1

test sparsersb
% ...
help sparsersb
% ...
Matrices created in this way support transparently nearly all of the operations supported by native GNU Octave sparse (see the GNU Octave's Sparse Matrices documentation).
If you are using a multi-core CPU and your matrices are significantly large, with sparsersb. you shall get a much faster matrix-vector multiplication than with Octave's sparse.
See the build instructions section on how to build sparsersb and the download section on where to get it.

A simple C Program using the Sparse BLAS Interface

Let us consider a simple program creating a sparse matrix and multiplying it by a vector.
Create a source code program file (say, mysbc.c) with a text editor: 
#include <rsb.h>         /* for rsb_lib_init */
#include <blas_sparse.h> /* Sparse BLAS on the top of librsb */
int main(const int argc, char * const argv[])
{
	blas_sparse_matrix matrix = blas_invalid_handle; /* handle for A */
	const int nnzA=4;	/* number of nonzeroes of matrix A */
	const int  nrA=3,ncA=3;	/* number of A's rows and columns */
	const int  incX=1,incB=1; /* spacing of X, B entries */
	int   IA[]={ 0, 1, 2, 2 }; /* row    indices*/
	int   JA[]={ 0, 1, 0, 2 }; /* column indices*/
	double VA[]={ 11.0, 22.0, 13.0, 33.0  }; /* coefficients */
	double X[]={  0.0,  0.0,  0.0 };
	double B[]={ -1.0, -2.0, -2.0 };
	rsb_err_t errval = RSB_ERR_NO_ERROR;
	if((errval = rsb_lib_init(RSB_NULL_INIT_OPTIONS)) != RSB_ERR_NO_ERROR) return -1;
	matrix = BLAS_duscr_begin(nrA,ncA);                             /* begin matrix creation */
	if(matrix==blas_invalid_handle) return -1;                    /* check */
	if( BLAS_ussp(matrix,blas_lower_symmetric) != 0 ) return -1;  /* set symmetry property */
	if( BLAS_duscr_insert_entries(matrix, nnzA, VA, IA, JA) != 0 ) return -1; /* insert some nonzeroes */
	if( BLAS_duscr_end(matrix) == blas_invalid_handle ) return -1; /* end creation */
	if( BLAS_dusmv(blas_no_trans,-1,matrix,B,incB,X,incX)) return -1;    /* X := X + (-1) * A * B  */
	if( BLAS_usds(matrix)) return -1;                              /* destroy matrix */
	if((errval = rsb_lib_exit(RSB_NULL_EXIT_OPTIONS)) != RSB_ERR_NO_ERROR) return -1;
	return 0;
}
To compile and build, execute in your shell: 
export PATH=$PATH:$HOME/local/librsb/bin/
gcc `librsb-config --I_opts --cflags` -c mysbc.c
gcc -o mysbc mysbc.o `librsb-config --static --ldflags --extra_libs`
./mysbc
See the reference documentation for more complete examples.

A simple modern C++ Program using the rsb.hpp Interface

The C++ <rsb.hpp> header introduced in librsb-1.3 allows C++ programs like mysbp.c: 
#include <rsb.hpp>
#include <vector>
#include <array>

int main() {
  rsb::RsbLib rsblib;
  const int nnzA { 7 }, nrhs { 2 };
  const int nrA { 6 }, ncA { 6 };
  const std::vector<int>    IA {0,1,2,3,4,5,1};
  const int                 JA [] = {0,1,2,3,4,5,0};
  const std::vector<double> VA {1,1,1,1,1,1,2}, B(nrhs * ncA,1);
  std::array<double,nrhs * nrA> C;
  const double alpha {2}, beta {1};

  // The three arrays IA, JA, VA form a COO (Coordinate) representation of a 6x6 matrix
  rsb::RsbMatrix<double> mtx(IA,JA,VA,nnzA); // Declarations of IA,JA,VA are all accepted via <span>

  mtx.spmm(RSB_TRANSPOSITION_N, alpha, nrhs, RSB_FLAG_WANT_ROW_MAJOR_ORDER, B, beta, C);
}
and compile and build, execute in your shell: 
export PATH=$PATH:$HOME/local/librsb/bin/
g++ -std=c++20 `librsb-config --I_opts --cxxflags` -c mysbp.cpp
g++ -o mysbp mysbp.o `librsb-config --static --ldflags --extra_libs`
./mysbp
See the reference documentation for more complete examples.

A simple Fortran Program using the Sparse BLAS Interface

The same program (say, mysbf.f90) in Fortran: 
      PROGRAM MAIN
      USE rsb             ! Fortran interface to C's rsb.h
      USE blas_sparse     ! Fortran interface librsb's Sparse BLAS
      IMPLICIT NONE
      INTEGER :: istat = 0             ! error status variable
      INTEGER :: A                     ! BLAS sparse matrix descriptor
      INTEGER :: transt=blas_no_trans  ! no transposition
      INTEGER :: incX=1, incB=1        ! X, B vectors increment
      REAL(KIND=8) :: alpha=3
      INTEGER :: nnzA=4, nrA=3, ncA=3  ! nonzeroes, rows, columns of matrix A
      INTEGER :: IA(4)=(/1, 2, 3, 3/)  ! row    indices
      INTEGER :: JA(4)=(/1, 2, 1, 3/)  ! column indices
      REAL(KIND=8) :: VA(4)=(/11.0, 22.0, 13.0, 33.0/) ! coefficients
      REAL(KIND=8) :: X(3)=(/   0,    0,    0/)
      REAL(KIND=8) :: B(3)=(/-1.0, -2.0, -2.0/)
      TYPE(C_PTR),PARAMETER :: EO = C_NULL_PTR ! 
      istat = rsb_lib_init(EO)                 ! librsb initialization
      IF(istat.NE.0)STOP                       ! check
      ! the 'd' in 'duscr_begin' specifies 'double precision' type
      CALL duscr_begin(nrA,ncA,A,istat)        ! begin matrix creation
      ! in the following calls the type is implicit due to the arguments
      CALL ussp(A,blas_lower_symmetric,istat)  ! set symmetry property
      CALL uscr_insert_entries(A,nnzA,VA,IA,JA,istat) ! insert some nonzeroes
      CALL uscr_end(A,istat)                   ! end matrix creation
      CALL usmv(transt,alpha,A,B,incB,X,incX,istat) ! X := X + (3) * A * B 
      ! please note the 'usds' is type independent
      CALL usds(A,istat)                       ! destroy matrix
      istat = rsb_lib_exit(EO)                 ! librsb finalization
      END PROGRAM MAIN
Note how default indexing in Fortran is 1-based.
To compile and build, execute in your shell: 
export PATH=$PATH:$HOME/local/librsb/bin/
gfortran `librsb-config --I_opts --fcflags` -c mysbf.f90
gfortran -o mysbf mysbf.o `librsb-config --static --ldflags --extra_libs`
./mysbf
The same program can be rewritten using the native RSB interface via the ISO C Binding feature: 
      PROGRAM MAIN
      USE rsb
      USE ISO_C_BINDING
      IMPLICIT NONE
      INTEGER,TARGET :: errval
      INTEGER :: transt = RSB_TRANSPOSITION_N  ! no transposition
      INTEGER :: incX = 1, incB = 1        ! X, B vectors increment
      REAL(KIND=8),TARGET :: alpha = 3, beta = 1
      INTEGER :: nnzA = 4, nrA = 3, ncA = 3     ! nonzeroes, rows, columns of matrix A
      INTEGER,TARGET :: IA(4) = (/1, 2, 3, 3/)  ! row    indices
      INTEGER,TARGET :: JA(4) = (/1, 2, 1, 3/)  ! column indices
      INTEGER :: typecode = RSB_NUMERICAL_TYPE_DOUBLE
      INTEGER :: flags = RSB_FLAG_DEFAULT_MATRIX_FLAGS+RSB_FLAG_SYMMETRIC
      REAL(KIND=8),TARGET :: VA(4) = (/11.0, 22.0, 13.0, 33.0/)     ! coefficients
      REAL(KIND=8),TARGET :: X(3) = (/   0,    0,    0/)
      REAL(KIND=8),TARGET :: B(3) = (/-1.0, -2.0, -2.0/)
      TYPE(C_PTR),TARGET  :: mtxAp=C_NULL_PTR
      TYPE(C_PTR)  :: mtxApp = C_NULL_PTR
      errval = rsb_lib_init(C_NULL_PTR)                             ! librsb initialization
      IF(errval.NE.0) STOP "error calling rsb_lib_init"
      mtxAp = rsb_mtx_alloc_from_coo_begin(nnzA,typecode,nrA,ncA,flags,C_LOC(errval)) ! begin matrix creation
      errval = rsb_mtx_set_vals(mtxAp,C_LOC(VA),IA,JA,nnzA,flags) ! insert some nonzeroes
      mtxApp = C_LOC(mtxAp)                                         ! create a matrix pointer pointer
      IF(errval.NE.0) STOP "error calling rsb_mtx_set_vals"
      errval = rsb_mtx_alloc_from_coo_end(mtxApp)                   ! end matrix creation
      IF(errval.NE.0) STOP "error calling rsb_mtx_alloc_from_coo_end"
      errval = rsb_spmv(transt,C_LOC(alpha),mtxAp,C_LOC(X),incX,C_LOC(beta),C_LOC(B),incB) ! X:=X+(3)*A*B 
      IF(errval.NE.0) STOP "error calling rsb_spmv"
      mtxAp = rsb_mtx_free(mtxAp)                                   ! destroy matrix
      IF(errval.NE.0) STOP "error calling rsb_mtx_free"
      errval = rsb_lib_exit(RSB_NULL_EXIT_OPTIONS)                  ! librsb finalization
      IF(errval.NE.0) STOP "error calling rsb_lib_exit"
      PRINT *, "rsb module Fortran test is ok"
      END PROGRAM MAIN
And linked the usual way. Just make sure you enabled the Fortran headers installation at configure time.
The advantage of using this interface is that it matches exactly the C version; this applies to the documentation as well.

Change Log, NEWS

In a nutshell: For a detailed list of user visible changes, see the NEWS.TXT file.
To get announcements, subscribe to this low-traffic Mailing List.

License & Co.

Mailing List

Mailing list (low traffic) is at: http://sourceforge.net/p/librsb/mailman/, mostly for librsb release announcements.

Contact the Author, Support


Links to Open Source Projects Using librsb

This webpage participates to the Viewable With Any Browser Campaign. Valid HTML 4.01 Transitional GNU Public License Frequently Asked Questions Partnership for Advanced Computing in Europe (PRACE)